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Abstract 

We consider simulations of a 2-dimensional gas of hard disks in a rectangular container and 
study the Lyapunov spectrum near the vanishing Lyapunov exponents. To this spectrum are 
associated "eigen-directions", called Lyapunov modes. We carefully analyze these modes 
and show how they are naturally associated with vector fields over the container We also 
show that the Lyapunov exponents, and the coupled dynamics of the modes (where it exists) 
follow linear laws, whose coefficients only depend on the density of the gas, but not on aspect 
ratio and very little on the boundary conditions. 

1 Introduction 

In this paper, we study the Lyapunov spectra of two-dimensional hard-disk systems and, in par- 
ticular, the associated "Lyapunov modes" El IS- Recently, this topic has received considerable 
attention |l3l|4l|5l|6ll3[8l, and a lot of progress has been made in the understanding of the issues 
involved. In the present work we synthesize and expand the results found earlier. In particular, 
we completely classify these modes and give a simple interpretation of their dynamics, in partic- 
ular for systems with arbitrary aspect ratio. We furthermore present new simulations in support 
of this classification. 

The Lyapunov exponents describe the rates of exponential growth, or decay, of infinitesimal 
phase-space perturbations, and are taken to be ordered according to A^^^ > A^^^ > ■ • ■ > A'^^\^ 

'i? is at most the dimension of the phase space. Note that A*-*' stands for different Lyapunov exponents, while we 
will use Aj when we consider them with multiplicity. 
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Because of the Hamiltonian nature of the problem, they come in conjugate pairs, 

As is well-known, is always a Lyapunov exponent for such systems and therefore, as a con- 
sequence, £ is odd. At any point ^ in phase space, the tangent space TX(Q decomposes into a 
sum 

TXiO = E^'KO © ■ ■ ■ © E''\0 , 

where E^^\^) is the (linear) space of those perturbations of the initial condition ^ whose growth 
rate is A*^^^ for the forward dynamics, and — A^-'^ for the time reversed dynamics. This decompo- 
sition is called the Oseledec splitting. 

We say that the Lyapunov exponent A^-'^ is d-fold degenerate if dim E^^\$^) = d. It should be 
noted that, when the Lyapunov exponents are rf-fold degenerate, only the subspace corresponding 
to all d of them is well defined^. The main idea of our full classification of Lyapunov exponents 
and their modes is based on this simple observation. 

The Lyapunov modes are defined as follows: at time t = 0, we take n orthogonal tangent 
vectors at ^ and, by applying to them the tangent-space dynamics^ for a long-enough time t, 
map them onto n vectors which, generally, are not orthogonal but still span an n-dimensional 
subspace S^{t). If, instead of n, we consider only n — 1 vectors, they similarly span an n — 1 
dimensional subspace S^-iit), such that S„_i{t) C 5'„(t). The Lyapunov mode is a unit vector 
in S„{t) which is orthogonal to the space S„^i{t). We will give a precise definition in the next 
section, explain the algorithmic aspects in Sect. 13. H and relate the modes to Oseledec's subspaces 
in Sect.lT2l 

The study of Lyapunov modes ^ IH has revealed interesting spatial structures which we 
will define later but which come in two types: localized structures associated with the large 
positive and negative Lyapunov exponents, and smooth delocalized structures of wave-like type 
for exponents close to zero |lj,8J. The exponents associated with the latter are degenerate and 
give rise to a step-like appearance of the Lyapunov spectrum as is shown, for example, in Fig.[Tl 

The discovery of these structures has led to several studies |l3l|4l|5l|6ll3 which go some way 
in explaining their origin and their dynamics. In this paper, we show how they are related to the 
symmetries of the container in which the particles move (including the boundary conditions). 
At the same time we obtain a classification of the degeneracies of the Lyapunov exponents near 
zero. This classification allows to view the so-called "mode dynamics" from a new geometrical 
perspective. 



^In this, and many other aspects, the theory of Lyapunov exponents is very similar to that of matrices. 
^See Section|2lfor details. 
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Figure 1: Lyapunov spectrum for = 780 hard disks at a density p = N/(L^Ly) = 0.8 in a 
rectangular periodic box with an aspect ratio Ly/L^ = 0.867. The insert provides a magnified 
view of the mode regime. I is the Lyapunov index numbering the exponents. 



The paper starts with a summary of results, passes through a precise definition of the Lya- 
punov modes, and then describes them as vector fields. In Sect.HJthese vector fields are classified, 
and their dynamics is studied in Sect. |5l The last sections deal with the density dependence and 
with possible hydrodynamic aspects. 

2 Notation and summary of results 

We consider a system of hard disks of diameter a and mass m moving in a two-dimensional 
rectangular container with sides of lengths and Ly. At this point, we do not need to specify 
the boundary conditions."^ The phase space of such a system is 

X = R2^x([0,LJx[0,L^])^, 

and a phase point ^ in X is 

^ = (P,q) = (Pi,---,Pn,(1i,---,(1n) ■ 

When necessary, we will write = {qi ^, Qi^y) to distinguish the two position components of 
particle i, and similarly for the momenta. The dynamics of the system is that of free flight, inter- 
rupted by elastic binary collisions. If is the state of the system at time 0, then = $*(^o) 
the state at time t, where $* : X ^ X defines the flow. 




"^We will consider reflecting and periodic boundary conditions. 
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Apart from issues of differentiability, which will be addressed when we describe the numer- 
ical implementation in Section |3l we call D$* the tangent flow. Informally speaking, it is a 
AN X AN matrix of partial derivatives and can be thought of as the first-order term (in e) in an 
expansion of the perturbed flow, 

+ = n^o) + ^ ^"^X ■ + 0{e') . 

The vector 5^ lies in the tangent space TX (at ,^q) to the manifold X; in our case TX(^q) = R^^. 
We next invoke the multiplicative ergodic theorem of Oseledec ||9l [T0|. To do so, we need 
ergodicity of the Liouville measure. Without further knowledge, we assume this for our case. 

Theorem: There exist an integer C. and numbers = such that for almost every 

^ G X (with respect to the Liouville measure) the tangent space splits into 

TX(0 = E^^Ko © E^^Ko © • ■ ■ © E^'Ko , 

with the property: 

lim ^ log II D<l>*L5f|| = ±\^j^ 

for 5^ e E'^^XO- The spaces E'^^XO are covariant: D$*|^ E^\0 = ^(^^($*(0)- 

The dimension of E^^\^) is called multiplicity of the exponent A'^-'^ In general the E^^^ are 
not orthogonal to each other. 

We also use the notation 

Ai > A2 > . . . > A, > . . . > 

to denote the Lyapunov exponents repeated with multiplicities, where the index is referred to as 
the Lyapunov index. This notation is more adequate for describing numerical methods of mea- 
surement, in which tangent-space dynamics is probed by a set of vectors. With these notations 
the relation between the A^^^ and the Aj is given by 

where /(^^ = S^^ + ■■■ + S'J\ Choosing j e {1, ...,£} and letting F^'J^ = E^^^ © ■ ■ ■ © E'-^\ we 
define the S^^ Lyapunov modes associated with A*^-'^ as any orthogonal spanning set of the space 

= (F(J"i))^ n F(^) . (1) 

A finer decomposition of this spanning set will be obtained when we describe the algorithm used, 
see Sect. 13. II The subspaces M^^^ are very similar to the subspaces E^^^: they have the same di- 
mension and also satisfy F^^^ = M^^^ © ■ ■ ■ © M^^\ However, they are not identical, because of 



Notation and summary of results 



5 



the orthogonality constraint in (HI). This will be explained in Sect. 13.21 

In this summary, we focus on rectangular boxes with periodic boundary conditions. Narrow 
systems {Ly < 2a) or systems with reflecting boundaries have a very similar Lyapunov spec- 
trum, but some exponents found in the periodic case are either absent, or appear with smaller 
multiplicities. We shall treat such systems in Sect. l4.3l but concentrate, until then, on the "gen- 
eral" periodic case. However, it should be noted that systems with reflecting boundaries give 
important information on the relation between the vanishing and the small Lyapunov exponents 
(see Ex.|6lbelow). 

The Lyapunov exponents near zero are found to be proportional to the wave numbers of the 
system. 



Remark: This result, as well as all results mentioned below, are to be understood in the limit of 
an infinite number of disks, at fixed density. In particular, we omit higher order terms (in k) in 
most of our statements. 

The Lyapunov exponents have the following properties: 

Lyapunov Spectrum: For gases of hard disks, the Lyapunov exponents near zero are fully de- 
termined by two (positive) constants, c-^ and c^^. For small-enough n = (n^, Uy), these exponents 
lie on two straight lines. ^ 

1 ) Transverse branch (T) .■ A = ic^ k^, with multiplicity 4 (2 if either or Uy is zero) 

2) Longitudinal branch (L) .• A = ±Cl k^, with multiplicity 8 (4 if either or Uy is zero) 

The multiplicities of both branches will be explained by simple geometric observations in Sect. 14.21 
While the linear laws resemble the square roots of the eigenvalues of a Laplacian in the box, we 
have no explanation beyond those already given in Ref. lISj. (The square root is related to the 
symplectic nature of the problem.) Additional degeneracies arise in square systems and may 
accidentally appear also for specific aspect ratios Ly/L^. 

We next explain how to visualize a mode EllHI • Fix ^ = (p, g) G X, and let 5^ be a Lyapunov 
mode of TX(Q. The vector 5C, = (Sp, 5q) has AN components, 2N associated with the momenta 
and 2N with the positions. Consider, for example, the q components, ^g^, . . . , 5qj^, where each 

^For both the longitudinal and transverse modes the linear k dependence of A is only the first term of an expansion 
in powers of A: |8|, A = c/c + Cjfc^ . . .. For positive A, is positive, but small. For a second order calculation in k. 




(2) 



see 1711. 
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Figure 2: Modes for the 780-disk system with periodic boundaries characterized in Fig.[T] Left: 
Transverse mode T(l,l) belonging to A1549. Right: Longitudinal mode for A1535 belonging to an 
LP pairLP(l,l). 



5qj is in (corresponding to the infinitesimal x and y displacements of g^, j = 1, . . . , A^). By 
drawing the perturbation vectors 5qj at the positions qj of the particles in the box, one obtains a 
field of vectors as is shown in Fig. El For dense-enough fluids we obtain a vector field in every 
point of the box by interpolating between the particles. 

It has been observed that for those Lyapunov exponents close to zero, this vector field is well 
approximated by trigonometric functions of the spatial coordinates x and y [l'.'W\. In particular, 
the number of nodes (n^, Uy) of the vector field determines a wave number k^. We then say that 
5^ is a mode of wave number k^. Our second main result is: 

Mode Classification: The subspaces M^\^) (defined in (|7])j belonging to Lyapunov exponents 
\^^^ close to zero fall into two categories: 

1 ) Transverse branch: modes associated with a Lyapunov exponent ±Crp k^ are divergence-free 
periodic fields of wave number 

2) Longitudinal branch: The modes associated with the Lyapunov exponent ±Cl k^ are of two 
types: 

( i) Half of them are irrotational periodic vector fields of wave number k^ ( called L-modes ); 

(ii) The others are scalar modulations with wave number k^ of the momentum field, (called 
P or momentum modes). 



The L and the P-modes turn out to be paired: the P-mode has components p^A{q^) and the corre- 
sponding L-mode has components V A{q^), where A is a scalar function, see ©. In the following, 
we refer to an L-mode and its corresponding P-mode as an LP pair. We will give more details 
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on this relation below. Results in this direction have been obtained in Refs. ||4l|6l|7|. 

An interesting question, which is related to some hydrodynamic aspects of the fluid, is an 
apparent propagation of the L and P-modes in physical space. It is a consequence of the motion 
of the tangent vectors in the subspaces M^^\^^), to which we refer as mode dynamics. It will be 
described in detail in Section |5l We note that this motion is not only determined by the tangent 
flow, but also by the re-orthonormalization process part of the algorithm for the simulation. Mode 
dynamics was observed early on ElIHl, but attempts to compute the propagation velocities are 
still scarce and, at present, only work for low densities |El|7l. Here, we give a more precise 
definition and provide numerical results. 

Mode Dynamics: In an M^-'^ space of Longitudinal and P-modes, the mode dynamics couples 
LP pairs. When restricted to the two-dimensional subspace spanned by a given LP pair, it reduces 
to a rotation at constant angular velocity which is proportional to the wave number k^, 

where v has the dimension of a velocity. 

In the remainder of the paper we state these results more precisely and give details about 
how they are obtained. They are of two types: First, a theoretical description of the modes, 
which is based on the symmetries of the system. Second, a detailed account of the numerical 
algorithms necessary to substantiate our claims (the difficulty being the decomposition of the o?^-' 
dimensional spaces M^^\Q). 

3 Tangent-space dynamics 

The dynamics of hard disks consists of phases of free flight interrupted by instantaneous elastic 
collisions. We denote the map for free flights of duration r by F"", and the collision map by C. 
Then, the evolution of an initial state, ^q, is given by 

= o • ■ ■ o C o o C o F^i ^0 , (3) 

where r^, . . . , r„ are the time intervals between successive collisions. The tangent-space dynam- 
ics of an infinitesimal perturbation 6^^ of is given by the tangent map of the flow ^ : 

where denotes the state just before the k^^ collision, and = C{^^) is the state immediately 
after.^ Here, DF'^L and DC L are AN x AN symplectic matrices. For the sake of simplicity, the 
flow is also called $*: namely $*(.Co) = ^i^d we often write 6^)- = D$*|^ ■ 6C,q. By convention, 
if t is a collision time, 5^^ denotes a tangent vector immediately before that collision. 

^We disregard problems of differentiability which appear for the (rare) tangent collisions. 
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3.1 Numerical procedure and Lyapunov modes 

Extensive numerical simulations are used to establish the classification and dynamics of the 
modes. Here, we briefly summarize our algorithm described in detail in Refs. lfT2t 21 and provide 
a precise definition of the Lyapunov modes we observe. It is well known that the exponential 
growth rate of a typical /c-dimensional volume element is given by + ■ ■ ■ + A^, 

Ai + • • • + A, = lim ^ log \\5^l A---ASm- (4) 

If the tangent vectors 5^^, . . . , are linearly independent at time zero, they will remain so, 
because the tangent flow is time reversible. Orthogonality is not preserved by the tangent flow, 
because, generally, D$* is not an orthogonal matrix. However, Q still holds if the vectors 
6^1, . . . , 5^^ are replaced by a (Gram-Schmidt) orthogonalized set of vectors 5r]j, . . . , Sr]^.^ It 
takes thus the simpler form 

1 ^ 

Ai + --- + Afc= lim-^log||Hl- (5) 

Since Q holds for every k < AN, the exponent A^, turns out to be equal to the growth rate of the 
/c*'^ orthogonalized vector 5ri^. 

The numerical method of Benettin et al. lfT3l and Shimada et al. lfT4ll - and indeed any al- 
gorithm - is based on this construction, although its basic objects are not the 57]^ vectors. As 
time increases, they all would get exponentially close to the most-unstable direction, become nu- 
merically indistinguishable, and diverge. Instead of orthogonalizing once at time t, the tangent 
dynamics is applied to a set of tangent vectors, which are periodically replaced by an orthog- 
onalized set, that we denote by ^7^^ . . . , . The modified dynamics is therefore that of an 
orthogonal frame lfT2l l2ll. One assumes that Q is still valid if 6r]^ is replaced by S'j^. The k^^ 
Lyapunov mode (or Lyapunov vector) at time t is, by definition, the vector S'j^, and it is associ- 
ated with the exponent A^. 

Our study starts with the observation that we consider a system with non-trivial Lyapunov 
exponents. To guarantee the convergence of the numerical algorithm, additional properties of the 
dynamical system are needed: Namely, that the system has well-defined local stable and unsta- 
ble subspaces associated with every Lyapunov exponent (close to zero). Results in this direction 
have been obtained for hard-disk systems in lfT5lfT6lfT7ll . A stronger property, hyperbolicity, has 
been recently proved for hard disk systems with randomly chosen masses in lfT8lfT9ll20ll . We as- 
sume here that these results hold for our system as well. As already advocated in Ref. fTO|, what 
matters from a physicist's point of view is that the numerical studies behave as if this were true. 
Under the above assumptions, the A;*'^ mode will align with the corresponding unstable subspace. 
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In other words, measured modes will be orthogonal spanning sets of the M^^^ subspaces defined 
in Eq. ©. 

The algorithm we use in our numerical work lfT2l l2ll is based on the principles just outlined. 
We restrict our considerations to hard-disk systems without external interaction. Since there is 
no potential energy, the dynamics is the same at any (total) energy, up to a rescaling of time. The 
natural unit of time of the system is {{ma'^N)/ K)^/"^ . Throughout, reduced units are used, for 
which the particle mass m, the disk diameter a, and the kinetic energy per particle, K/N, are 
unity. The density, defined by p = N/V and the aspect ratio, defined hy A = Lyj L^, are the only 
relevant macroscopic parameters. Here, V = L^Ly is the area of the (rectangular) simulation box 
whose sides are and Ly in the x and y directions, respectively. All our numerical examples 
are for densities p < 0.8 characteristic of dense or dilute (if p < 0.1) hard disk gases. The results 
are insensitive to the time between successive Gram-Schmidt re-orthonormalization steps. 

3.2 Modes and Oseledec subspaces 

We have remarked in Sect.|2lthat the spaces E''^^ of Oseledec' theorem are in general not identical 
to the spaces of the modes M^^^ defined in the Eq. ([T]). Here, we explain this difference. 

The covariant subspaces E^^^^ . . . , E^^^ of the Oseledec splitting are obtained as follows ifTOll : 
the multiplicative ergodic theorem (for reversible systems) states that the matrices 

1 

A±(0 = lim (D<I>±*|J D<I>±*L)2|<I 

exist with probability one. The eigenvalues of are exp(A*^^) > . . . > exp(A^^^) and the 
eigenvalues of A_ are exp(— A*^^) > . . . > exp(— A^^^). Since both A^ and A_ are symmetric, 
their eigenspaces define two orthogonal decompositions of the tangent space 

where U'£^ is the eigenspace of A_|_ associated to exp(±A'^-'-'). Note that the subspaces are 
pairwise orthogonal but in general not covariant. However, for j E {!,...,£}, the subspaces 

U^f®...®Uf and t/i^^ © . . . © f/^'^ 

are covariant. They are, respectively, the subspace of the i — j + 1 most stable directions of 
A_,_ and the subspace of the j most unstable directions of A_ (or equivalently its j most stable 
directions in the past). For every ^, the invariant subspaces E^^^ are then given by 

E^'^ = (f/i^) © ... © u^l^) n (t/f © ... © uf) . 

The -E^^ spaces are covariant but in general not orthogonal. One can show that 

?7i^^ © . . . © U^I^ = E(^) © ... © = F^^') , 
where F''^^ was introduced in Sect.|2l Using the definition dl]), one easily verifies that M ^•'^ = U'^\ 
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Remark 1 In our numerical method, the modes are obtained at the phase points $*(.Co) 
large t, i.e., at the end of integrated trajectories. It is therefore consistent that they should contain 
information about the past and not the future, that is about A and not A^. 

In general, U'f \ f/+^ and E''^'' are different. Hard-disk systems have a special property: the 
three spaces coincide for the Lyapunov exponent A*^*"^ = 0, where m = |(£ + 1) is the middle 
index. We shall denote by J\f the covariant subspace J\f = E^"^\ also called the null subspace. In 
Sect. 14.11 we give the explicit form of J\f. Using the explicit Jacobians for the particle collisions 
and for the free- streaming motion, as given in llfT2l l22l l2l. one can check that the orthogonal 
complement to the null subspace, J\f^, is also covariant. From this statement, one can show that 
J\f = f/<™) = U^_^\ Our measured modes M^™^ = [/i™^ therefore span the null space. Since this 
argument seems to fail for j ^ m, we can only conclude that modes of M^^^ are perturbations 
whose exponential growth rate is at least 

3.3 Two simplifications 

Two additional properties greatly simplify the classification of the modes: the symplecticity of 
the tangent flow and an observed property is proportional to Sp". The symplecticity of the 
tangent flow means that 



D<I>*|JJD<I>*L = J, (6) 



where the 4N x 4N matrix J is defined as 

J 



Sq) \1 OjUqJ { 6p 

As a consequence of Q, both A^ and A_ are symplectic (see e.g., EJ) and, as is well known, 
Lyapunov exponents of Hamiltonian systems come in pairs A^-'^ A*^^"-'"'"^^ = — A'^-'^ of equal mul- 
tiplicity. The matrix J relates the U^'' subspaces by 

U2\0 = JUt'^\0 . (7) 

Therefore, it is sufficient to measure the Lyapunov exponents and the modes for the positive part 
of the spectrum. Note that J is not the (derivative of the) velocity reversal {5p, 5q) ^— {—6p, 6q) 
and does not change the sign of time. 



The simulations of our system exhibit an additional structure which considerably simplifies 
the analysis of Lyapunov modes near zero. Unstable perturbations (A > 0) have the (approxi- 
mate) form 

5p = C({, X)5q, C > , 
while the stable (A < 0) perturbations are of the form 



-C'i^A)Sp = Sq , C'>0. 
'With the reduced units defined above, positions and velocities are dimensionless. 
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Here, C and C are numbers. In the left panel of Fig.|3]we demonstrate that the perturbations 5q 
and 5p associated with each Lyapunov mode are nearly parallel or anti-parallel for large A^. The 
Equation dV]) also suggest C = C, which is well verified by the simulations. 




-I I 1 ^ I ^ 

1000 2000 3000 

/ 



Figure 3: Value of cos(0) = (6q ■ Sp)/(\6q\ ■ \6p\), as a function of the Lyapunov index / 
for an instantaneous configuration of the system characterized in Fig. [B Here, is the angle 
between the 2A^-dimensional vectors of the perturbation components of all particle positions and 
velocities for a Lyapunov vector For the small positive exponents, for which / < 2N — 2 = 
1558, this angle vanishes, for the small negative exponents, for which / > 2N + 3 = 1563, it 
is equal to vr. For the six zero modes, 1558 < / < 1563, the angle lies between these limiting 
values. 

We can therefore restrict our classification to the 6q part of the modes corresponding to 
positive exponents. If we can associate every measured exponent A to a given 5q, then we will 
know that the exponents A and —A correspond to the modes {5p, Sq) with Sp = CSq and Sp = 
—C^^6q, respectively. 

3.4 Tangent vectors as vector fields 

The components of a tangent vector 5^ = {5p, 5q) are the perturbation components of the po- 
sitions and velocities of all particles. As a graphical representation of such a vector, we show 
in the left panel of Fig. |4|the instantaneous positional perturbations of all particles at their po- 
sitions in physical space, where the arrows indicate the directions and strengths. It belongs to 
a Lyapunov exponent A 1540 indicated by the enlarged circle in Fig. |5l A qualitatively identical 
figure is obtained, if the positional displacements of all particles are replaced by their momen- 
tum displacements, as explained in Sect. 13.31 Thus, this figure is a complete representation of 
the Lyapunov vector 5E, = {5p, 5q) belonging to Xi^^^ in Fig.|5l The transverse modes for Xi^^^ 
on the left panel of Fig. |2l and for A 1543 at the bottom left panel of Fig. |6l are other examples of 
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Figure 4: The transverse mode T(l,l) for of Fig.|5l Left: Interpretation as a vector field; 
Right: Alternative representation as periodic spatial patterns of the position perturbations Sx^ and 
of the particles, which emphasizes the wave vector parallel to a diagonal of the simulation 
box. 



transverse modes belonging to the same degenerate exponent. 

We interpret the left panel of Fig. |4| as a two-dimensional vector field which - up to a 
constant phase - is well described by 

V (Py(x, y)J V -"2 sin(A;^ x) cosiky y) ) ' 

where k^ = and ky = and ai, are two constants. For this reason, we assign the node 
numbers (n^., n^) = (1, 1) to this mode. 

To be more precise, let r = (x, y') E [0, L^) x [0, Ly). We say that a two-dimensional smooth 
vector field (p = (ip^, ipy) over the position space, is sampled by the infinitesimal displacements 
5q of the disks at their reference positions q, when 

V?(gP = 6qj , for all j = 1, . . . , . 

Thinking of a tangent vector in terms of an associated vector field is meaningful if there are "suf- 
ficiently" many particles to sample the field on a typical length scale of its variation. Once this 
condition breaks down for larger exponents, as it will for large-enough k, the modes disappear^, 
and so do the steps in the Lyapunov spectrum. In the following we use a notation which does not 
distinguish between tangent vectors and their two-dimensional vector fields. 

''and the sampled vector field is ill-defined. 
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Figure 5: Lyapunov spectrum and "dispersion relations" for the 780-disk system characterized in 
Fig.IU and the density p = 0.8. Left: Lyapunov exponents are ordered by size and repeated with 
multiplicities. The specially-marked point corresponds to the transverse mode shown in Fig.|4| 
Right: Lyapunov exponents as a function of their wave number. The respective labels L and T 
refer to the longitudinal and transverse branches. 



4 Observation and description of the modes 

In this section, we describe the Lyapunov spectrum near and the corresponding modes, as they 
are measured in numerical experiments. In the two following subsections, we will explain how 
these modes can be understood on the basis of symmetry breaking of "zero modes". 

We illustrate our assertions with the system already introduced in Fig. [H which contains 
= 780 particles in a rectangular periodic box with an aspect ratio L^/L^ = 0.867. It corre- 
sponds to a hard-disk gas with a density p = 0.8, slightly below the fluid-to-solid phase transition 
density llT2l . The left panel of Fig.|5]provides a magnified view of the smallest positive Lyapunov 
exponents for this system.^" The exponents, ordered by size and repeated with their multiplic- 
ities, are plotted as a function of their index. Degenerate exponents with a multiplicity d > 2 
appear therefore as "steps". To account for the wave-like appearance of the modes, we associate 
with each Lyapunov vector a wave number k, as in Q, where the non-negative integers 

n = (n^, Uy) count the nodes in the respective directions.'^ 

When the small Lyapunov exponents are plotted as a function of their corresponding wave 
number, they all lie on two curves, sometimes referred to as "dispersion relations" This is 

'"The conjugate negative exponents are not shown, see ||2||8]. 

"in many cases, the number of nodes is easy to determine, as in Fig.|4] To facilitate an objective identification of 
the wave vectors modes associated with larger exponents, which are noisier, Fourier-transforms are used. 
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demonstrated on the right panel of Fig.|5l On this plot, degenerate exponents are represented by 
a single point. For reasons discussed below, the upper branch is called longitudinal (L), and the 
lower transverse (T). It is experimentally found that for a given wave number, the multiplicity of 
the L branch is twice that of the T branch, as mentioned already in Sect. El 

4.1 Vanishing Lyapunov exponents 

We start our description with the six modes associated with the six vanishing Lyapunov expo- 
nents, commonly referred to as zero modes. Four of them are induced by the homogeneity of 
space, and two are consequences of the homogeneity of time. They span a six-dimensional sub- 
space J\f{0 of the tangent space, TX(Q, at any phase point These 6 zero modes play a 
fundamental role in understanding the nature of the modes associated with Lyapunov exponents 
close to zero. 

In order to show which symmetries give rise to the zero modes, we list in Table [l] the six 
corresponding elementary transformations. This defines the six zero modes 5^^ to in a nota- 
tion that separates the x and y components of 5p and 5q. The vectors and 5^2 correspond to a 
perturbation of the total momentum in the x and y directions, 5^^ and 5^^^ to an (infinitesimal) uni- 
form translation of the origin, 5^^ to a change of energy, and to a change of the origin of time. 



Tr ansf or mat ion 


Generator 


(P, q) ^ ^x^Py + ^^Ax^Qy) 
(P, ^ {Px,Py,(lx + ^^^(ly) 

(p, q) ^ (Px,Py,qx,qy + ^'^) 

(p, q) ^ (Px + ^Px, Py + ^Py, qx, qy) 

(P, q) ^ (Px, Py,qx + ^Px, qy + ^Py) 


5^1 = (1,0, 0,0) 
5^2 = (0,1, 0,0) 
5^3 = (0,0, 1,0) 
5^4 = (0,0, 0,1) 

^^5 = (Px,P?;, 0,0) 

= (0,0,p^.,py) 



Table 1: Central subspace for the vanishing Lyapunov exponents. Notation : 1 = (1, 1, . . . , 1), 
and = (0, 0, . . . , 0). All vectors have AN components. 

One can explicitly check that the subspace Span{5^^, . . . , 5^g} is covariant and that its vec- 
tors have a sub-exponential growth (or decay). It thus coincides with the null space A/'(0 of 

Sect. Ha 

Remark 2 The space A/'(0 can be further decomposed into three covariant subspaces = 
Span{(5^i, (^^s), A/"^ = Span{5^2) '^^4} -^p = Span{(5^5, 5^^}, each of which independently 
satisfies the properties listed for A/'(0- In systems with reflecting boundaries, only A/^ is present, 

'^It is important to keep in mind that M{0 really depends on ^, see below. 
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while and My are absent, because they are related to translation invariance. If only the x 
direction is periodic |6|, the space of zero modes is reduced to A4 © M^. 

4.2 Longitudinal, Transverse and P-modes 

Following our argument of Sect. 13.31 we need only describe the 5q part of the modes. Therefore, 
we consider the three transformations (see Table [l]) 

■ (QxJ, Qyj) ^ (QxJ, QyJ + ^) (8) 

■ (QxJ, Qyj) ^ (QxJ + £Px,j, QyJ + ^Py,j) , 

associated to zero modes, and claim the following : 

Modes Classification I: Modes of wave number are scalar modulations of (0) with wave 
number k^. More precisely, they are obtained by replacing e with e A(q^ ,j, q^j) in (0), where the 
real scalar function A is of the form 

A(x, y)= ^i,m exp k^x + m k^y)) . 

\£\=n^,\m\=ny 

The space of such modulations has dimension 4 in general and dimension 2 if either or Uy 
vanishes. 

Examples We consider n = (1,0), (0,1) and (1,1). We use the notation = cos(k^x), 
= cos(k^x), Cy = cos(kyy) and Sy = sm(kyy). A basis of the space of functions with wave 
number k^ is shown in Table |2l Of course, this choice fixes a constant phase for the sine and 
cosine functions. 



n 


Function A 


dim 


(1,0) 


^xy ^x 


2 


(0,1) 


Cy, Sy 


2 


(1,1) 


^x^yi ^x^y, C^Sy, S^Sy 


4 



Table 2: Functions with wave number k^ 



If the Lyapunov exponent were a function of k^ only, we would expect 12-fold degeneracy 
in Fig. 121 (resp. 6-fold if either or Uy = 0), since each of the three perturbations of (EJ) can be 
modulated by the four (resp. two) functions of Table|2l However, this degeneracy is broken into 
an 8-1-4 (resp. 4-1-2) structure: 

Mode Classification II: The Lyapunov vectors of wave number k^ have two possible Lyapunov 
exponents: |A| = or |A| = Cl jA^nl. corresponding to the Transverse or Longitudinal 

branch ofFig.\5\ The modes for these two branches are obtained as follows: 
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1 ) Transverse branch: transverse modes are obtained by combining the modulations of 6^^ 
and (5^4 in a divergence-free vector field. We denote by T(n) the space of such vector fields. 

2) Longitudinal branch: 

( i) Longitudinal modes are irrotational vector fields one obtains by combining the modula- 
tions of 6^^ and 6^^. We denote the corresponding space by L(n). 

(ii) P-modes are modulations of5^Q, and we denote the corresponding subspace by P(n). 

The three subspaces T(n), L(n) and P(n) have dimension 4 (or dimension 2 if either or Uy 
vanishes). We denote by LP(n) = L(n) © P(n) the subspace corresponding to the Longitudinal 
branch. It has dimension 8 (4 if either or n vanishes). 



Remark 4 The divergence and curl of a vector field if = ((p^, cpy) are, of course, 

V ■ = d^<f^ + dy<fy , V A = d^fy - dyf^ . 

Since every two-dimensional vector field uniquely decomposes into a sum of a divergence-free 
and an irrotational vector field, the spaces L(n) and T(n) span all possible two-dimensional 
vector fields of wave number A;„. 

There is a simple way to build P(n), L(n) and T(n) from the scalar modulations of Table El 
If A is a modulation of wave number k^, then we have 

pA eP(n), VAGL(n), V A A G T(n) , (9) 

where by definition 

This construction is also useful because it naturally defines what we shall call LP pairs, by 
which we denote a field of L(n) and a field of P(n) originating from the same scalar modulation, 
as in Q. We show below that LP pairs play an important role. Indeed, when only some of the 
modes are present because of the boundary conditions, LP pairs are never broken: both fields are 
present, or both are missing (see Sect. l4.3l) . We shall also see in Sect.|5lthat the dynamics of the 
modes mostly takes place between LP pairs. 

Example 5 For the three modes of lowest wave number. Table |3l lists a basis of T(n), L(n) and 
P(n). Fields are given in a non-normalized form to keep notation short. Corresponding L(n) 
fields and P(n) fields are LP pairs. Fig.|6lprovides examples for T and L-modes. 

Modes of L(n) and T(n) are wavelike perturbations of the position space, and, therefore, are 
similar to the modes that appear in hydrodynamics (see Sect.|7l). In particular, when either or 
Uy vanish, the fields of L(n) and T(n) are, respectively, longitudinal and transverse to the wave 
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n 



Basis of T(n) 



Basis of L(n) 



Basis of P(n) 



(1,0) 
(0,1) 

(1,1) 



o\ /o 



/' V 



Cx^y J V S.j.Sy 



oy' V 
o\ /o 



^y 



y 



k '^x^yX ( k ^x^y 

1 '1 

k ^x^^y' Vfc ^x^y 



ky ^x^y \ / ky ^x^^y 

~T^x(^y^ ^~T^x^y 



PyJ \Py 



PyJ ^ \Py, 



Px\ fPx 
Py \Py) ''"''^ 



Px\ (Px 

\Py) """"^ 



Table 3: Decomposition of n-modes for rectangular systems with periodic boundaries. 



vector. This observation is the reason for the names of the two branches in Fig.|5l To stay in line 
with this now-accepted terminology, we keep it also for the case ■ Uy ^ 0. P-modes are more 
complex than the other modes, because they depend not only on the positions of the perturbed 
particles but also on their velocities. 

At this point we are able to characterize the spectrum and its multiplicities with only two 
constants, Cl and c^, which have the physical dimension of velocities. We shall see in Sect. |6l 
that these velocities depend on the density of the system, but are insensitive to the system size, 
the boundary type, and the aspect ratio. Therefore, it is tempting to think of Cl and as thermo- 
dynamic velocities. However, as discussed in Sect.|6l no obvious interpretation could be found 
so far. 

4.3 Other aspect ratios and boundary conditions 

We summarize here some observations which concern different boundary conditions and degen- 
eracies. 

i) The Lyapunov spectrum is more degenerate for square systems, for which n)-^{n n y 
In this case the multiplicities are doubled with respect to the general case, for which ^ 
Uy. Other "accidental" degeneracies may occur: for instance, parameters can be found for 
which Ct |fc(i,i)| = Cl |A;(i o)I- 

ii) Systems with reflecting boundaries Q develop only a subset of the modes encountered so 
far, which can be found by the following simple, and obvious, rules: 

• The fundamental wave vectors are = and ky = j- (not 2n). Here, and Ly 
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Figure 6: Snapshots of Lyapunov modes for the periodic 780-disk system of Sect. IH Left, 
from top to bottom: Transverse modes T(1,0), T(0,1), and T(l,l) belonging to A1557, A1555, and 
Ai54g, respectively, of Fig. |5l Right, from top to bottom: Vector fields for longitudinal modes, 
which belong to the LP pairs LP(1,0), LP(0,1), and LP(1,1), and which are associated with the 
respective exponents A1553, A1545 and A1535, of Fig.|5l 
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are the effective box sizes, obtained from the actual side lengths by subtracting one 
particle diameter a. 

• The fields of T(n) and L(n) have to satisfy Dirichlet conditions, namely to be tangent 
to the boundary: 

V5x(0, y) = ^x(Lx, y) = and (py(x, 0) = (py(x, Ly) = . 

If expressed in terms of sines and cosines, this means that may contain sm(k^x) but 
not cos{k^x), and so on. 

• As for periodic boundary conditions, when an L-mode is present, its paired P-mode is 
always present. 

iii) Hybrid systems with a reflecting boundaries in one direction and periodic boundaries along 
the other behave as expected |l6|: the fundamental wave vectors are chosen according to the 
boundary type, and the Dirichlet conditions are only applied to one component of the field. 

iv) Narrow systems, for instance those with a < Ly < 2a and ^ a,^^ only develop modes 
with Uy = ll23l . This follows, since a vector field varying along the y axis cannot be 
sampled by a single particle. Therefore, the Lyapunov spectrum of such a system is greatly 
simplified, since the modes are restricted to L{n^, 0) and T(n^,, 0). 

Example 6 Table |4| shows which modes of Table |3] satisfy the Dirichlet condition and are thus 
present in a system with reflecting boundaries. We stress that such a system has only two 
vanishing Lyapunov exponents which are associated with S^^ and of Table [U Therefore, 
modes appear which are not modulations of zero modes of the system. The crucial observation 
here is that even if one of the fundamental symmetries is broken by the boundary condition, 
the modulation, as defined in Eq. Q, may still satisfy that boundary condition. For example, 
if we have reflecting boundaries on the walls {x = 0} and {x = L^}, then any perturbation 
A(x, y) = siB.{mk^x)B(y) will be acceptable, whereas, of course, A(x, y) = cos{mk^x)B{y) 
would not. 

Systems with reflecting boundaries and narrow systems are easier to study numerically, be- 
cause the multiplicities of the L, P and T spaces are smaller than in the periodic rectangular case. 
In particular, the LP(1,1) space has only dimension 2 when the boundaries are reflecting. For 
that reason, we illustrate some of the issues below also with the "reflecting- wall version" of the 
780-disk system introduced in Fig.lU 



Recall that a is the diameter of the disks. 
''^For such a choice of basis vectors the origin of the coordinate system is at the bottom-left corner of the simula- 
tion box. 
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n 


T(n) 


L(n) 


P(n) 


(1,0) 


none 


Co) 




(0,1) 


none 


C) 




(1,1) 






ill) ''"''^ 



Table 4: Mode decomposition for rectangular systems with reflecting boundaries. 



4.4 How to measure P-modes 

Contrary to L(n) and T(n), the tangent subspace P(n) depends on the state that it perturbs. 
Moreover, modes of P(n) are not really vector fields, because the velocities of the particles in a 
typical configuration do not depend smoothly on position. In order to "see" a P-mode, we face 
two problems: 

i) A typical measured vector of an LP(n) space is a superposition of vectors of L(n) and P(n). 

ii) Even when a mode of P(n) is isolated, it is not smooth and does not "look like" a vector 
field over the box. 

We explain our solution to the first problem with the simplest possible example: the LP(1, 0) 
space for a rectangular system with reflecting boundaries (see Sect. l4.3l) . This space has dimen- 
sion two and is defined by the two normalized spanning vectors 

1 (SA ^ (Px\ „ , ^ 



where and Z2 are normalization constants. At any time t, we have two measured modes, and 
ip'^ (vectors with 2N components, since only the 6q part is observed), whose span is (numerically 
very close to) LP(1, 0). Therefore, there are constants a, b, c, d with 

if^ = '^^a + ^^b, ^^^^ 
Lf^ = c + d . 

Since all these vectors are normalized, one should also have d^ + b"^ = + d'^ = 1. However, 
Eq. (fTTT) is over-determined: four constants have to satisfy AN equations. Numerically, we use 
a least-square method to find the best values for a,b,c and d, which we denote by a, (3, 7, 5P 

'^It is also possible to use a simple projection a — ip^ ■ -0^, etc. This method assumes that and Lp^ are 
orthogonal vectors, which they only are approximately. 
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Figure 7: Example for the P-mode reconstruction for the 4-dimensional LP(1,0) space of a system 
with periodic boundaries. Only half of the components are shown, namely those corresponding 
to s^. The components are similar. Top row: x and y components of the reconstructed P-mode 
^ ps^. Bottom row: the fields /p^ and (^y /py are wavelike again. Note that in the top row 
one can recognize the sinusoidal envelope. 



Thus, the measured modes, ij)^ and ij)'^, are decomposed according to 

(p^ = iIj^ a + ip"^ (5 , 

where the vectors and are the best-possible LP-pair, tp^ and (p^, reconstructed from ex- 
perimental data. 

Now we can deal with our second problem: whereas the vector is easily recognized as a 
vector field, the vector is not. However, from (fTOb we deduce that 
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or, more precisely, 



Therefore, both ^p^/p, 



X 




(12) 



Remarks 



i) This procedure readily extends to higher-dimensional LP(n) spaces, for instance to the 
four-dimensional LP(1, 0) space of a periodic system. We take this case as an example to 
illustrate in Fig. |7] our method of mode reconstruction: from the four measured modes the 
analogue to Eq. (fTTTl generates four spanning vectors, namely two L-modes (one of which 
is shown at the top-right position of Fig. |4|) and two P-modes. The x and y components 
of one of the P-modes, namely (^^ ^ ps^, are shown in the top row of Fig. IT] No smooth 
functions are recognized. However, after dividing by the momentum components, p^, ^ and 
Py j, as required by (fT2b . the figures for the fields (p^/p^ and 0^/Py in the bottom row of 
Fig. clearly display the expected -dependence. 

ii) Division by p^ j or Py^^ in (fT2b is numerically unstable when particle j has a very small 
momentum along a coordinate axis. We avoid this by multiplying instead by p^j/iplj + £), 
where e <^ 1. One could also ignore particle j in this case and only sample the field at those 
points where both p^ j and Pyj are not too small. 

5 Dynamics of the modes 

Remark The interested reader can look up animated pictures on the web at the address 
|http : //theory . physics . unige . ch/modes^. 

In this section we turn to the "dynamics of the modes" and study what has been called the 
velocity of the longitudinal modes |l2l|5l|7|. This might clarify hydrodynamic theories El El, 
which are mostly based on a partial classification of the modes. In numerical simulations, the 
transverse modes are stationary in space and time: although the particles move, the vector-field 
of the mode does not. In other words, at any instant of time, the vector field of a T-mode does 
not move (up to a small jitter due to numerical noise). 

For longitudinal modes, however, one seems to observe ll24ll a propagation in the direction of 
the wave vector. Using the geometrical picture developed above, we can interpret this motion 
and also explain why no propagation is observed for the LP dynamics in systems with reflecting 
boundaries as demonstrated in Fig.[8]below. 



when either or n vanish, for instance in the L(l, 0) or L(0, 1) space. 
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Figure 8: LP(1,1) dynamics for a system of 780 disks in a box with an aspect ratio 0.867, a 
density 0.8, and with reflecting boundaries. Left: coordinates of the measured fields in the 
"standard basis" of LP(1,1). The (unit) circle is nearly reached, showing that indeed ip'^ and ip'^ 
span the same subspace as (p^ and (p^ . Center and right: the measured fields ip^it) and ip'^(t). 
The rows from top to bottom are consecutive snapshots separated by time steps of At = 3.20, 
which corresponds to a phase shift of 7r/2 in the LP(1,1) rotation. 
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Since multiplicities are not so essential here, we illustrate the interpretation for the (simpler) 
LP(1, 0) space of a rectangular system with reflecting boundaries (dimension 2). In Sect. 14.41 we 
pointed out that at any given time t, where the state is the measured modes ip^it) and ip'^it) 
are combinations of the two spanning vectors (^^(^t) and (^^(^t). We re-write Eq. (fTTTl in matrix 
form, but now with explicit time dependence: 



Therefore, the dynamics of the modes reduces to that of the 2x2 matrix Q(t). In our simulations 
we keep the spanning vectors orthonormal. Since (in the experiment) the two fields ^p^ and ^p^ 
are also (nearly) orthogonal'^, the matrix Q{t) is close to a rotation matrix (that is c ~ —b, 
d a, and + 6^ ~ 1). Therefore, the dynamics in the two-dimensional subspace is well 
described by a phase 0(t) = arctan(6(t)/a(t)). 

LP-dynamics in 2 dimensions: The matrix Q{t) is a rotation with constant angular velocity 
uj^. This velocity is proportional to the wave number k^, namely 



Here, v is the product of a frequency and a wavelength and, therefore, is a velocity. According 
to the simulations, v depends only on the density of the system. 

Example 7 In Fig.|8]we demonstrate this rotation in the LP(1,0) space of a system with reflecting 
boundaries. Snapshots of the two measured fields at consecutive times show the rotation of the 
two vectors between the L and the P direction. 

A similar rotation has been found in narrow systems in |l6| and explained in |l5l |7| in the low 
density limit using a Boltzmann-equation approach. 

Remark 8 The velocity v = uj^/l^^nl has been interpreted earlier ll24l as the phase velocity of 
a traveling wave in physical space. In Appendix |B] we demonstrate how the two interpretations 
can be reconciled. The definitions given above allow us to apply the same concepts also to the 
LP-dynamics of systems with reflecting boundaries, although they do not show traveling but 
standing waves. 



'^The scalar product ip^ ■ tp^ = Il^^-^cos(k^qj ^)sm(k,^q^ ^)pj ^ a priori does not vanish. However, as the 
simulations show, it is of the same order as ^f^i cosik^q^ ^) siii(fc^gj which is also small and non-vanishing due 
to the uneven spacing of the particles. 




a(t) 
c{t) 




0(t) = uj^t = vk^t . 
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Figure 9: Rotations of various LP pairs for the 780-disk system of Fig. [B Left: Time dependence 
of the phase. For clarity, the phases for different modes are separated by multiples of2'K/\k^\. 
Right: The fluctuations around the constant velocity (in percent of 2tx). For details we refer to 
the main text. 



Next, we consider the dynamics for the general case of a 2(i-fold degenerate LP(n) space. 
As explained in Sect. IH an LP(n) space is defined by a spanning set of d longitudinal modes, 
Lp\, . . . , and d P-modes, y^f , . . . , y^^, where each pair {Lp\, Lp^) is an LP pair. The following 
description is valid for any type of boundary conditions: 

Dynamics in LP(n): The dynamics in LP(n), when restricted to a two-dimensional subspace 
spanned by an LP pair 

Span{v7^,(^n, 

is a two-dimensional rotation at a constant angular velocity icu^, where uo^ = v k^, with v in- 
dependent ofn. 

Example 9 The LP dynamics for higher-dimensional spaces is illustrated in Fig. |9l In the fol- 
lowing we concentrate onLP(l,l). This space has dimension 8, that is four LP pairs. Generically, 
a measured mode has a non-vanishing projection onto all four LP pairs, and four different phases 
can be defined. The four time series of phases for LP(1,1) in Fig.|9lbelong to different projections 
of the same mode. 

6 Influence of geometry and system size 

In this section we study what influence the density, aspect ratio, and boundary conditions have 
on Cx, Cl, and v. The density dependence is significant. Unfortunately, we do not have any 
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Figure 10: Slopes of the transverse and longitudinal branches and Cl of Fig.|5l and of phase 
velocity v. The simulations are for a system containing N = 780 particles in a rectangular 
periodic box with a fixed aspect ratio Ly/L^ = 0.867. Left: c-p and Cl as a function of the 
particle density p. The smooth lines are polynomial fits added to guide the eyes. Although 
the fits have no theoretical basis, we provide the fit parameters for convenience: Cp = 0.790 + 
0.970p + 0.785p2 + 6.24p^ and = 0.902 + l.OSOp - O.OSSp^ + 3.85p\ Right: Almost linear 
relations between all three quantities. 



explanation for this fact. In particular, comparisons with the sound velocity ll25l . with the mean 
free path, and with similar quantities, do not suggest simple relationships. Thus, these questions 
have to await further studies. 

The results of our simulations are summarized in Fig.[lOl In the left panel, Cp and Crp, defined 
in Sect. |2l are shown as functions of p, and on the right panel v and c-p are plotted as functions 
of Cl. There are two observations which are of interest: First, as seen in the left panel, Cp and Cp 
cross for lower densities. Thus, the intuitively natural conjecture, Cp > Cp, is not supported by 
the numerical evidence. Second, the three "velocities" Cp, Cp, and v are almost linearly related. 
In particular, the phase velocity v agrees rather well with Cl for fluid densities p < 0.6 (see the 
full line in the right panel of Fig. [TOb . thus lending support to referring to the curves A(fc) as 
"dispersion relations". 

On the left panel of Fig.[TT]we plot the longitudinal and transverse dispersion relations for 
a moderately-dense gas with a density p = 0.4. The full and open points are for systems with 
periodic and reflecting boundaries, respectively. In all cases, A was determined from the lowest 
step of the Lyapunov spectrum for the system, for which the aspect ratio varied between 0.5 and 
1, and the particle number between 400 and 800. The figure demonstrates that periodic and re- 
flecting boundary conditions give the same A(A;). We have experimentally verified (but not shown 
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Figure 11: Simulation results for a hard-disk gas with a density p = 0.4. The full and open 
points refer to periodic and reflecting boundary conditions, respectively. Left: The smallest 
positive Lyapunov exponent for T and L-modes, respectively, as a function of the wave-number 
k. The curves for periodic and reflecting boundary conditions agree. To the second order in k, 
a fit to the data gives = 1.422A; + 0.55P and = 1.412A; + O.OSQP. Right: The slopes, 
and Crp, are plotted as a function of the aspect ratio A of the simulation box. For a given A we 
compare points with the same k to eliminate the influence of the nonlinearity of the dispersion 
relations. This means that for periodic boundaries the linear extensions of the simulation box, 
and Ly, are twice those of the respective reflecting box. The data for periodic and reflecting 
boundaries agree to within numerical accuracy. 



here) that the same is true also for larger and lower densities. 

The situation becomes a little more complicated when we consider the dependence of the 
slopes, Cl and Crp, on the aspect ratio A, as we do in the right panel of Fig.[TT] The dispersion 
curves are not strictly linear in k, as we have pointed out already in a footnote in Sect. |2l A fit 
to the points in the left panel of Fig. [TT] reveals that the term proportional to k"^ is much larger 
for L than for T. To eliminate this nonlinearity in a comparison of the slopes for reflecting and 
periodic systems with a given A, modes with the same values for k are used in the right 
panel of Fig.[TT] With this precaution the figure demonstrates that the slopes of the dispersion 
relations do not depend on the boundary conditions in any significant way, as long as the box 
does not degenerate to a narrow channel El. Except for this particular case, the nonlinearity 
of the dispersion curves has no noticeable influence on the classification and description of the 
Lyapunov modes in this paper and has been accordingly ignored. 

'^The smallest wave number k for periodic boundaries is given by 2tt/L, where L is a dimension of the box. For 
reflecting boundaries, k is given by n/L, where L is an effective box size, namely the size reduced by a particle 
diameter a. 
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7 Hydrodynamic equivalent of the modes 

There is a general expectation that the Lyapunov modes should be related to the hydrodynamic 
behavior of the system, and several papers point in this direction |l3l@l|5l|6l|. It should be noted, 
however, that none of these studies has reached a totally convincing interpretation, and, further- 
more, it is obvious from this body of work that the LP-modes are more difficult to explain than 
the T- modes. Here, we add to this a simple calculation which might be helpful in the future: we 
determine how the modes perturb the hydrodynamic fields or, in other words, what would be the 
hydrodynamic equivalent of the modes we measure. The results, given in Table |5l have quite a 
simple form, but do not reproduce the usual modes of hydrodynamics. 

Consider a general transformation T of the one-particle phase space /i = [0, L^) x [0, Ly) x 
given by 

^ f r I— > r' = r + e 5^{r,v) 
' y V ^ v' = V + e 5r](r, v) ' 

and / a probability density over /i. If e is infinitesimal, then the new probability density is 
f = f + e6f, where 

Sf = -f (V, ■ 5^ + V, ■ 5v) - VJ ■ 5^ - V J ■ 5r] , (13) 

see Appendix lAl Since we study equilibrium dynamics, we assume, furthermore, that / is the 
Boltzmann distribution'^. By Sect. 13. 31 we also have 5r] = C5^, and, therefore, (fT3t simplifies to 

= -f (V, ■ 6^ + CV, -S^-Cv 60 ■ 

We define the variations of the three hydrodynamic fields (density p, momentum u, and energy 
E) by 

6p(r)= / dv 5f(r,v) , 5u(r) = / dv 5f(r,v)v , SE(r)= / dv 5f(r,v)\v\'^ . 



mode 


5^ 


bp 


bu 


bE 


T 


V AA 





V A A 





L 


WA 




VA 




P 


V A 





VA 


A 



Table 5: Modulations of the hydrodynamic fields (multiplicative constants are omitted) 

For each of the three types of modes, one can compute these quantities as a function of the 
scalar modulation A introduced in Q. The results are given in Table |5l Since AA = —k\A, we 
note that the scalar fields bp and bE are proportional to the initial scalar modulation A. Note also 
that the energy field is only affected by the L and P-modes but not by the T-mode. 

'^with /crT = 1 as in the simulations 
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8 Conclusions 

The picture developed in this paper is for two-dimensional hard disks, but the method is suffi- 
ciently geometric and general to allow easy extensions to other systems: 

For example, the generalization to three-dimensional hard disks is straightforward. The exis- 
tence of L and T-modes for this case has been confirmed by computer simulation ll25ll . Recently, 
Lyapunov modes were also found for two-dimensional soft-particle systems interacting either 
with a Weeks-Chandler- Anderson potential [i26J or a Lennard- Jones potential ll27l |2^ I . It would 
also be interesting to extend the work to, say, circular geometries. Another extension concerns 
linear molecules such as hard dumbbells in a periodic box dllll- In this case two qualitatively 
different degrees of freedom play a role, translation and rotation. The existence of modes has 
already been demonstrated in this case. 

We have ended our wanderings through the rich landscape of Lyapunov modes. To summa- 
rize, we have carefully identified and analyzed the modes, giving a beginning of a theoretical 
classification. Furthermore, we have seen that the Lyapunov exponents and the phase velocity 
of the LP-modes seem to be functions of the density alone. In particular, they are practically in- 
dependent of the aspect ratio of the box (and, where applicable, also insensitive to the boundary 
conditions). 
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A Transformation of the one-particle distribution 

Let / be a given distribution. For F C /i, we have 



The transported probability F' is defined by F'{V') = F(T), where F' = T(T). We shall compute 
its density /'. In the integral 



F(T) = Prob(one particle G F) 



dr dv f(r, v) . 



r 




dr dv f(r, v) 



(14) 



we change the variables to (r', v') = T(r, v). To first order in e, we have 




r' t-H> r = r' — e 6^(r', v') 
ti' h- > y = — £ 6ri{r', v') 
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so (1141) becomes 

F(D = j^^ dr' dv' det [ BT-\,^^, ] fir' - e S^{r\ v'), v' - e 5r]{r\ v')) . (15) 
The determinant is 

det[Dr-|,J = det(I^-^y^ -^Y. ) 

V -eT),d7] Id - £:D„(5r7 y ^/y (16) 

To first order in e, we obtain with (fTSb and (fT6b 

F(r) = F(rO - £ dr' dv' fir', v') [V, -S^ + V,- 5t]\^,^^, 
- e j^^ dr' dv' [VJ ■ 5i + VJ ■ ^r^],, , 

which is equivalent to (flSt . 

B Traveling waves of LP dynamics 

We consider the LP(1, 0) space of a rectangular system with periodic boundary conditions, 
defined by the four spanning vectors 

We take an initial tangent vector 

^^0 = «<^s^in + &<^L ■ (17) 

After a time r = 271 /{Auj^, the vector is transformed to 

Assume that a and h are more or less equal. Then (flTt will resemble a sinus, with much "noise" 
due to the P component, while dTSt will look more like a cosine. In the dynamics leading from 
(flTt to dTSt a kind of "traveling wave" is therefore visible, which seems to cover a distance 
27r|/Cn|^^ in a time 2nuj^^, thus "moving " at velocity v = c^n/l^nl- actual simulations, 
we cannot expect typical vectors to have a phase difference of | between their cp]^^^ and (^^os 
components, as in our example. Therefore, the observed wave displacement as seen in ll25l l24ll 
has the shape of "steps" in a space-time diagram, with an average slope equal to v. 
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